c2.prop=function(data,phati,mi,Mi,n,N){
  data=as.name(data); N=N; phati=phati
  mi=mi; Mi=Mi; qhati=1-phati
  n=n; fpc=(1-n/N)
  phat=sum(Mi*phati)/sum(Mi)
  Mbar=sum(Mi)/n
  s2r=sum((Mi^2*(phati-phat)^2))/(n-1)
  vhat=fpc*(1/(n*Mbar^2))*s2r+(1/(n*N*Mbar^2))*sum(Mi^2*(1-mi/Mi)*(phati*qhati/(mi-1)))
  bound=2*sqrt(vhat)
  lower=phat-bound; upper=phat+bound
  lowerpct=sprintf("%.2f%%",lower*100); upperpct=sprintf("%.2f%%",upper*100)
  cat("","\n","Results from 2-stage Cluster sample: Data =",data,"\n","\n","N =",N,"n =",n,"FPC =",fpc,"\n",
      "phat =",phat,"\n","Vhat =",vhat,"\n","Bound =",bound,"\n","Lower Bound =",lower,"Upper Bound =",
      upper,"\n","Lower Bound (%) =",lowerpct,"Upper Bound (%) =",upperpct,"\n","","\n")
  results=list(phat=phat,data=data,n=n,N=N,fpc=fpc,vhat=vhat,
               bound=bound,lower=lower,upper=upper,lowerpct=lowerpct,upperpct=upperpct,s2r=s2r)
}

# to use the function with its call:
# c2.prop(data,phati,mi,Mi,N)

# data: name of dataset, in quotes
# phati: vector (data) of proportions from clusters
# mi: vector (data) of elements from sample
# Mi: vector (data) of elements from population
# N: population size